--- title: Image Preprocessing keywords: fastai sidebar: home_sidebar nb_path: "nbs/02_preprocessing.ipynb" ---
{% raw %}
{% endraw %}

Especially for MRI the pixel values can vary between scanner types. This will lead to a very differnt range of pixel intensities across images in the same batch and the algorithms will have a hard time converging. Various techniques exist for rescaling exists of which some are implemented here.

{% raw %} 
@patch
def show_L(l:L, **kwargs): 
    "shows all Tensors in L"
    show_images_3d(torch.stack(tuple(l)), **kwargs)
    
@patch
def apply(l:L, fun):
    "applies a function to each element of L"
    return L(fun(item) for item in l)
{% endraw %} {% raw %} 
fns = Path('../data/series').ls()
images = L(TensorDicom3D.create(fn) for fn in fns)
images.show_L(nrow = 15)
{% endraw %}

Size correction

If working with data from differnt scanners of even different sequences, pixel spacing can differ. E.g. T2w images are usually in a higher resolution as DWI images, thus the size of each individual pixel is larger in DWI images than in T2W. Normalizing pixel size to a standartized value can help achiving better performance. However, only the H x W Dimension is rescaled, not the D Dimension.

{% raw %} 
from faimed3d.basics import TensorDicom3D, TensorMask3D # for compatibility with show_docs
{% endraw %} {% raw %}
{% endraw %} {% raw %}

TensorDicom3D.size_correction[source]

TensorDicom3D.size_correction(im:TensorMask3D'>), new_spacing=1)

{% endraw %} {% raw %}

TensorMask3D.size_correction[source]

TensorMask3D.size_correction(im:TensorMask3D'>), new_spacing=1)

{% endraw %} {% raw %} 
[im.get_spacing() for im in images]

[(0.6835939884185791, 0.6835939884185791, 1.0),
 (0.6835939884185791, 0.6835939884185791, 6.0),
 (0.6835939884185791, 0.6835939884185791, 2.0),
 (0.6835939884185791, 0.6835939884185791, 6.0),
 (0.6835939884185791, 0.6835939884185791, 6.0)]
{% endraw %}

Pixel rescaling

rescale intercept: 0028|1052
rescale slope: 0028|1053

{% raw %} 
from faimed3d.basics import TensorDicom3D # for compatibility with show_docs
{% endraw %} {% raw %}
{% endraw %} {% raw %}

TensorDicom3D.rescale_pixeldata[source]

TensorDicom3D.rescale_pixeldata(t:TensorDicom3D)

{% endraw %} {% raw %} 
im = images[0]
print(im.mean())
# the example MRI images have no rescale slope or intercept
im.set_metadata('0028|1052', 5.)
im.set_metadata('0028|1053', 0.5)
im = im.rescale_pixeldata()
im.mean()

TensorDicom3D(153.2178)

TensorDicom3D(81.6089)
{% endraw %}

Pixel Intensity Normalization

Mean / Max / Median scaling

{% raw %}
{% endraw %} {% raw %}

TensorDicom3D.mean_scale[source]

TensorDicom3D.mean_scale(t:TensorDicom3D)

Scales pixels by subtracting the mean and dividing by std. 0 pixels are ignored

{% endraw %} {% raw %}

class MeanScale[source]

MeanScale(p=1.0, **kwargs) :: RandTransform

A transform that before_call its state at each __call__

{% endraw %} {% raw %} 
images.apply(MeanScale()).show_L(nrow = 15)
{% endraw %} {% raw %}
{% endraw %} {% raw %}

TensorDicom3D.median_scale[source]

TensorDicom3D.median_scale(t:TensorDicom3D)

Scales pixels by subtracting the median and dividing by the IQR. 0 pixels are ignored

{% endraw %} {% raw %}

class MedianScale[source]

MedianScale(p=1.0, **kwargs) :: RandTransform

A transform that before_call its state at each __call__

{% endraw %} {% raw %}
{% endraw %} {% raw %}

TensorDicom3D.max_scale[source]

TensorDicom3D.max_scale(t:TensorDicom3D)

{% endraw %} {% raw %}

class MaxScale[source]

MaxScale(p=1.0, **kwargs) :: RandTransform

A transform that before_call its state at each __call__

{% endraw %}

Histogram scaling

However, just scaling wiht one mean (or median) and std (or IQR) might not be the optimal solution. Histgram scaling might be used for a more evendistribution of pixel values. Differnt functions are provided for histogram scaling:

Method adapted form Jeremy Howard

See the excelent Kaggle kernel form Jeremy Howard for an explanation

{% raw %} 
from torch import Tensor # for compatibility with show_docs
{% endraw %} {% raw %}
{% endraw %} {% raw %}

TensorDicom3D.freqhist_bins[source]

TensorDicom3D.freqhist_bins(t:Tensor'>), n_bins=100)

A function to split the range of pixel values into groups, such that each group has around the same number of pixels. taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

{% endraw %} {% raw %}

Tensor.freqhist_bins[source]

Tensor.freqhist_bins(t:Tensor'>), n_bins=100)

A function to split the range of pixel values into groups, such that each group has around the same number of pixels. taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

{% endraw %} {% raw %}

TensorDicom3D.hist_scaled[source]

TensorDicom3D.hist_scaled(t:Tensor'>), brks=None)

Scales a tensor using freqhist_bins to values between 0 and 1 taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

{% endraw %} {% raw %}

Tensor.hist_scaled[source]

Tensor.hist_scaled(t:Tensor'>), brks=None)

Scales a tensor using freqhist_bins to values between 0 and 1 taken from https://github.com/fastai/fastai/blob/master/fastai/medical/imaging.py#L78

{% endraw %} {% raw %}

TensorDicom3D.hist_scaled_pt[source]

TensorDicom3D.hist_scaled_pt(t:Tensor'>), brks=None)

same as fastai fucntion for PILDicom

{% endraw %} {% raw %}

Tensor.hist_scaled_pt[source]

Tensor.hist_scaled_pt(t:Tensor'>), brks=None)

same as fastai fucntion for PILDicom

{% endraw %}

The sample data from radiopaedia is already windowed an scaled, so the results presented here are not very impressive. However, if applied to MRI data, changes will be stronger.

{% raw %} 
def f(x): return x.hist_scaled()
images.apply(f).show_L(nrow = 15)
{% endraw %}

Comparison to SimpleITK Adaptive Histogram Equalization

The method from Jeremy Howard is more simple that methods provided in SimpleITK, however it is much faster.

{% raw %} 
im1 = TensorDicom3D(images[0][10:20])
{% endraw %} {% raw %} 
image = sitk.Cast(im1.as_sitk(), sitk.sitkFloat32)
output = sitk.AdaptiveHistogramEqualization(image, radius=[20]*3, alpha=0.3, beta=0.3) # ~ 30-60 seconds
TensorDicom3D(sitk.GetArrayFromImage(output)).show() 
im1.hist_scaled().show()
{% endraw %}

N4 Bias Field Correction

from official SimpleITK docs

The N4 bias field correction algorithm is a popular method for correcting low frequency intensity non-uniformity present in MRI image data known as a bias or gain field. The method has also been successfully applied as flat-field correction in microscopy data. This method assumes a simple parametric model and does not require tissue classification.

{% raw %} 
corrector = sitk.N4BiasFieldCorrectionImageFilter()
corrector.SetMaximumNumberOfIterations([3]*3)
output = corrector.Execute(sitk.Cast(im1.as_sitk(), sitk.sitkFloat32)) # needs float
TensorDicom3D(sitk.GetArrayFromImage(output)).show() #
{% endraw %}

faimed3d provides a wrapper class to apply both N4 Bias Field Correction and Adaptive Histogram Equalization, which basically just simplifies reading the data.

{% raw %}
{% endraw %} {% raw %}

class ImageCorrectionWrapper[source]

ImageCorrectionWrapper(n4_max_num_it=3, hist_radius=[5, 5, 5], hist_alpha=0.3, hist_beta=0.3, do_n4=True, do_hist=True, verbose=True)

{% endraw %}

Piecewise Linear Histogram Matching

[1] N. Laszlo G and J. K. Udupa, “On Standardizing the MR Image Intensity Scale,” Magn. Reson. Med., vol. 42, pp. 1072–1081, 1999.

[2] M. Shah, Y. Xiao, N. Subbanna, S. Francis, D. L. Arnold, D. L. Collins, and T. Arbel, “Evaluating intensity normalization on MRIs of human brain with multiple sclerosis,” Med. Image Anal., vol. 15, no. 2, pp. 267–282, 2011.

Implementation adapted from: https://github.com/jcreinhold/intensity-normalization, ported to pytorch (no use of numpy, works in cuda).

In contrast to hist_scaled, the piecewise linear histogram matching need pre-specified values for new scale and landmarks. It should be used to normalize a whole dataset.

{% raw %}
{% endraw %} {% raw %}

get_percentile[source]

get_percentile(t:Tensor, q:float)

Return the q-th percentile of the flattened input tensor's data.

CAUTION:

  • Needs PyTorch >= 1.1.0, as torch.kthvalue() is used.
  • Values are not interpolated, which corresponds to numpy.percentile(..., interpolation="nearest").

:param t: Input tensor. :param q: Percentile to compute, which must be between 0 and 100 inclusive. :return: Resulting value (float).

This function is twice as fast as torch.quantile and has no size limitations

{% endraw %} {% raw %}
{% endraw %} {% raw %}

get_landmarks[source]

get_landmarks(t:Tensor, percentiles:Tensor)

Returns the input's landmarks.

:param t (torch.Tensor): Input tensor. :param percentiles (torch.Tensor): Peraentiles to calculate landmarks for. :return: Resulting landmarks (torch.tensor).

{% endraw %} {% raw %}
{% endraw %} {% raw %}

find_standard_scale[source]

find_standard_scale(inputs, i_min=1, i_max=99, i_s_min=1, i_s_max=100, l_percentile=10, u_percentile=90, step=10)

determine the standard scale for the set of images Args: inputs (list or L): set of TensorDicom3D objects which are to be normalized i_min (float): minimum percentile to consider in the images i_max (float): maximum percentile to consider in the images i_s_min (float): minimum percentile on the standard scale i_s_max (float): maximum percentile on the standard scale l_percentile (int): middle percentile lower bound (e.g., for deciles 10) u_percentile (int): middle percentile upper bound (e.g., for deciles 90) step (int): step for middle percentiles (e.g., for deciles 10) Returns: standard_scale (np.ndarray): average landmark intensity for images percs (np.ndarray): array of all percentiles used

{% endraw %} {% raw %}
{% endraw %} {% raw %}

Tensor.piecewise_hist[source]

Tensor.piecewise_hist(image:Tensor, landmark_percs, standard_scale)

Do the Nyul and Udupa histogram normalization routine with a given set of learned landmarks

Args: input_image (TensorDicom3D): image on which to find landmarks landmark_percs (torch.tensor): corresponding landmark points of standard scale standard_scale (torch.tensor): landmarks on the standard scale Returns: normalized (TensorDicom3D): normalized image

{% endraw %} {% raw %}

class PiecewiseHistScaling[source]

PiecewiseHistScaling(landmark_percs=None, standard_scale=None, final_scale=None, p=1.0, slicewise=False, slice_dim=None, **kwargs) :: RandTransform

Applies theNyul and Udupa histogram nomalization and rescales the pixel values.

Args: input_image (TensorDicom3D): image on which to find landmarks landmark_percs (torch.tensor): corresponding landmark points of standard scale final_scale (function): final rescaling of values, if None is provided values are scaled to a mean of 0 and a std of 1. slicewise (bool): if the scaling should be applied to each slice individually. Slower but leads to more homogeneous images

Returns: If input is TensorMask3D returns input unchanged If input is TensorDicom3D returns normalized and rescaled Tensor

{% endraw %} {% raw %} 
standard_scale, percs = find_standard_scale(images)
standard_scale, percs

(tensor([313.1477, 343.0235, 348.7227, 353.8698, 359.4490, 364.1022, 367.5029,
         370.4589, 373.5394, 377.9434, 385.5532]),
 tensor([ 1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 99]))
{% endraw %} {% raw %} 
images.apply(PiecewiseHistScaling(landmark_percs=percs, 
                                  standard_scale=standard_scale)).show_L(nrow = 15)
{% endraw %} {% raw %}
{% endraw %} {% raw %}

standard_scale_from_filelist[source]

standard_scale_from_filelist(fns:Series'>))

calculates standard scale from images in a filelist

{% endraw %} {% raw %} 
from fastai.data.core import DataLoaders # for compatibility with show_docs
{% endraw %} {% raw %}
{% endraw %} {% raw %}

DataLoaders.standard_scale_from_dls[source]

DataLoaders.standard_scale_from_dls(dls:DataLoaders)

calculates standard scale from images in a dataloader

{% endraw %}